{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018 parts of this notebook are from [this Jupyter notebook](https://nbviewer.jupyter.org/github/krischer/seismo_live/blob/master/notebooks/Computational%20Seismology/The Finite-Difference Method/fd_ac1d.ipynb) by Heiner Igel ([@heinerigel](https://github.com/heinerigel)), Lion Krischer ([@krischer](https://github.com/krischer)) and Taufiqurrahman ([@git-taufiqurrahman](https://github.com/git-taufiqurrahman)) which is a supplemenatry material to the book [Computational Seismology: A Practical Introduction](http://www.computational-seismology.org/),  additional modifications by D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href=\"https://fonts.googleapis.com/css?family=Merriweather:300,300i,400,400i,700,700i,900,900i\" rel='stylesheet' >\n",
       "<link href=\"https://fonts.googleapis.com/css?family=Source+Sans+Pro:300,300i,400,400i,700,700i\" rel='stylesheet' >\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' >\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width */\n",
       "    width: 800px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 1000px;\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "margin-top:0.5em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256); \n",
       "    border-radius: 0px; \n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 110%;\n",
       "    width:680px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-style:regular;\n",
       "    font-weight: bold;    \n",
       "    font-size: 250%;\n",
       "    line-height: 100%;\n",
       "    color: #004065;\n",
       "    margin-bottom: 1em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: bold; \n",
       "    font-size: 180%;\n",
       "    line-height: 100%;\n",
       "    color: #0096d6;\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "\tfont-size: 150%;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: #008367;\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: 300; \n",
       "    font-size: 100%;\n",
       "    line-height: 120%;\n",
       "    text-align: left;\n",
       "    width:500px;\n",
       "    margin-top: 1em;\n",
       "    margin-bottom: 2em;\n",
       "    margin-left: 80pt;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    font-weight: regular;\n",
       "    font-size: 130%;\n",
       "    color: #e31937;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .5em;\n",
       "    margin-top: 1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'Source Code Pro', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "\t\n",
       "    \n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"], \n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FD stability and dispersion\n",
    "\n",
    "In the [last lesson](http://nbviewer.jupyter.org/github/daniel-koehn/Theory-of-seismic-waves-II/blob/master/03_Intro_finite_differences/2_fd_ac1d.ipynb) we developed a 1D acoustic FD modelling code. For the given modelling parameters, the code worked flawlessly and delivered modelled seismograms, which are in good agreement with the analytical solution. In this lesson we want to investigate how to choose optimum time steps dt and spatial grid point distances dx, to get stable and accurate FD modelling results. We start, by revisiting a simplified version of our 1D acoustic FD modelling code ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Import Libraries \n",
    "# ----------------\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from pylab import rcParams\n",
    "\n",
    "# Ignore Warning Messages\n",
    "# -----------------------\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [],
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Definition of modelling parameters\n",
    "# ----------------------------------\n",
    "xmax = 500 # maximum spatial extension of the 1D model (m)\n",
    "dx   = 0.5 # grid point distance in x-direction\n",
    "\n",
    "tmax = 1.001   # maximum recording time of the seismogram (s)\n",
    "dt   = 0.0010  # time step\n",
    "\n",
    "vp0  = 333.   # P-wave speed in medium (m/s)\n",
    "\n",
    "# acquisition geometry\n",
    "xr = 365.0 # receiver position (m)\n",
    "xsrc = 249.5 # source position (m)\n",
    "\n",
    "f0   = 25. # dominant frequency of the source (Hz)\n",
    "t0   = 4. / f0 # source time shift (s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison of numerical with analytical solution\n",
    "\n",
    "In the function below we solve the homogeneous 1D acoustic wave equation by the 3-point spatial/temporal difference operator and compare the numerical results with the analytical solution. To play a little bit more with the modelling parameters, I restricted the input parameters to dt and dx. The number of spatial grid points and time steps, as well as the discrete source and receiver positions are estimated within this function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     42
    ]
   },
   "outputs": [],
   "source": [
    "# 1D Wave Propagation (Finite Difference Solution) \n",
    "# ------------------------------------------------\n",
    "def FD_1D_acoustic(dt,dx):\n",
    "        \n",
    "    nx = (int)(xmax/dx) # number of grid points in x-direction\n",
    "    print('nx = ',nx)\n",
    "        \n",
    "    nt = (int)(tmax/dt) # maximum number of time steps            \n",
    "    print('nt = ',nt)\n",
    "    \n",
    "    ir = (int)(xr/dx)      # receiver location in grid in x-direction    \n",
    "    isrc = (int)(xsrc/dx)  # source location in grid in x-direction\n",
    "\n",
    "    # Source time function (Gaussian)\n",
    "    # -------------------------------\n",
    "    src  = np.zeros(nt + 1)\n",
    "    time = np.linspace(0 * dt, nt * dt, nt)\n",
    "\n",
    "    # 1st derivative of a Gaussian\n",
    "    src  = -2. * (time - t0) * (f0 ** 2) * (np.exp(- (f0 ** 2) * (time - t0) ** 2))\n",
    "\n",
    "    # Analytical solution\n",
    "    # -------------------\n",
    "    G    = time * 0.\n",
    "\n",
    "    # Initialize coordinates\n",
    "    # ----------------------\n",
    "    x    = np.arange(nx)\n",
    "    x    = x * dx       # coordinate in x-direction\n",
    "\n",
    "    for it in range(nt): # Calculate Green's function (Heaviside function)\n",
    "        if (time[it] - np.abs(x[ir] - x[isrc]) / vp0) >= 0:\n",
    "            G[it] = 1. / (2 * vp0)\n",
    "    Gc   = np.convolve(G, src * dt)\n",
    "    Gc   = Gc[0:nt]\n",
    "    lim  = Gc.max() # get limit value from the maximum amplitude\n",
    "    \n",
    "    # Initialize empty pressure arrays\n",
    "    # --------------------------------\n",
    "    p    = np.zeros(nx) # p at time n (now)\n",
    "    pold = np.zeros(nx) # p at time n-1 (past)\n",
    "    pnew = np.zeros(nx) # p at time n+1 (present)\n",
    "    d2px = np.zeros(nx) # 2nd space derivative of p\n",
    "\n",
    "    # Initialize model (assume homogeneous model)\n",
    "    # -------------------------------------------\n",
    "    vp    = np.zeros(nx)\n",
    "    vp    = vp + vp0       # initialize wave velocity in model\n",
    "\n",
    "    # Initialize empty seismogram\n",
    "    # ---------------------------\n",
    "    seis = np.zeros(nt) \n",
    "    \n",
    "    # Calculate Partial Derivatives\n",
    "    # -----------------------------\n",
    "    for it in range(nt):\n",
    "    \n",
    "        # FD approximation of spatial derivative by 3 point operator\n",
    "        for i in range(1, nx - 1):\n",
    "            d2px[i] = (p[i + 1] - 2 * p[i] + p[i - 1]) / dx ** 2\n",
    "\n",
    "        # Time Extrapolation\n",
    "        # ------------------\n",
    "        pnew = 2 * p - pold + vp ** 2 * dt ** 2 * d2px\n",
    "\n",
    "        # Add Source Term at isrc\n",
    "        # -----------------------\n",
    "        # Absolute pressure w.r.t analytical solution\n",
    "        pnew[isrc] = pnew[isrc] + src[it] / dx * dt ** 2\n",
    "                \n",
    "        # Remap Time Levels\n",
    "        # -----------------\n",
    "        pold, p = p, pnew\n",
    "    \n",
    "        # Output of Seismogram\n",
    "        # -----------------\n",
    "        seis[it] = p[ir]   \n",
    "        \n",
    "    # Compare FD Seismogram with analytical solution\n",
    "    # ---------------------------------------------- \n",
    "    # Define figure size\n",
    "    rcParams['figure.figsize'] = 12, 5\n",
    "    plt.plot(time, seis, 'b-',lw=3,label=\"FD solution\") # plot FD seismogram\n",
    "    Analy_seis = plt.plot(time,Gc,'r--',lw=3,label=\"Analytical solution\") # plot analytical solution\n",
    "    plt.xlim(time[0], time[-1])\n",
    "    plt.ylim(-lim, lim)\n",
    "    plt.title('Seismogram')\n",
    "    plt.xlabel('Time (s)')\n",
    "    plt.ylabel('Amplitude')\n",
    "    plt.legend()\n",
    "    plt.grid()\n",
    "    plt.show()            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dx   = 0.5 # grid point distance in x-direction\n",
    "dt   = 0.0010  # time step\n",
    "FD_1D_acoustic(dt,dx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the same result, we achieved in the last lesson. Now, you might get the smart idea to save some computation time by increasing the timestep dt. Let's try it ... "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dx   = 0.5 # grid point distance in x-direction\n",
    "#dt   = 0.0010  # old time step\n",
    "dt   = 0.0015023  # time step\n",
    "FD_1D_acoustic(dt,dx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Oops, maybe this idea was not so smart at all, because the modelling becomes **unstable**. Instead of increasing the time step dt, we could try to increase the spatial discretization dx to save computation time ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dx   = 0.5 # old grid point distance in x-direction\n",
    "dx   = 7.0 # new grid point distance in x-direction\n",
    "dt   = 0.0010  # time step\n",
    "FD_1D_acoustic(dt,dx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hmm, the accurracy of the FD modelling result compared to the analytical solution is clearly deterioated, when the spatial grid point $dx$ is increased. And why does the P-body wave becomes dispersive? More generally, how do I choose $dx$ and $dt$ without using a trial-and-error approach, which requires a lot of computation time, especially when considering 3D modelling. To understand the underlying problems, we will investigate the stability and numerical dispersion of the FD method in the next two sections."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stability of 1D acoustic wave equation finite difference approximation\n",
    "\n",
    "To analyse the stability of the finite difference approximation of the 1D acoustic wave equation:\n",
    "\n",
    "\\begin{equation}\n",
    " \\frac{p_{j}^{n+1} - 2 p_{j}^n + p_{j}^{n-1}}{\\mathrm{d}t^2} \\ = \\ vp_{j}^2 \\frac{p_{j+1}^{n} - 2 p_{j}^n + p_{j-1}^{n}}{\\mathrm{d}x^2},\n",
    "\\end{equation}\n",
    "\n",
    "we use an approach introduced by the famous mathematician and pioneer of computational sciences [John von Neumann](https://en.wikipedia.org/wiki/John_von_Neumann). For the **von Neumann Analysis**, we assume harmonic plane wave solutions for the pressure wavefield like:\n",
    "\n",
    "\\begin{equation}\n",
    "p = exp(i(kx-\\omega t)),\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "with $i^2=-1$, the wavenumber $k$ and circular frequency $\\omega$. Using the discrete \n",
    "spatial coordinates:\n",
    "\n",
    "$x_j = j dx,$\n",
    "\n",
    "and times \n",
    "\n",
    "$t_n = n dt.$\n",
    "\n",
    "We can calculate discrete plane wave solutions at the discrete locations and times in eq. (1), for example at grid point j and time n:  \n",
    "\n",
    "\\begin{equation}\n",
    "p_j^n = exp(i(kjdx-\\omega n dt)),\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "or at grid point j and time n+1:\n",
    "\n",
    "\\begin{align}\n",
    "p_j^{n+1} &= exp(i(kjdx-\\omega (n+1) dt))\\nonumber\\\\\n",
    "&= exp(-i\\omega dt)\\; exp(i(kjdx-\\omega n dt))\\nonumber\\\\\n",
    "&= exp(-i\\omega dt)\\; p_j^n,\\nonumber\\\\\n",
    "\\end{align}\n",
    "\n",
    "or at the grid point j and time n-1:\n",
    "\n",
    "\\begin{align}\n",
    "p_j^{n-1} &= exp(i(kjdx-\\omega (n-1) dt))\\nonumber\\\\\n",
    "&= exp(i\\omega dt)\\; exp(i(kjdx-\\omega n dt))\\nonumber\\\\\n",
    "&= exp(i\\omega dt)\\; p_j^n.\\nonumber\\\\\n",
    "\\end{align}\n",
    "\n",
    "Similar approximations can be estimated for time n at the spatial grid points j+1:\n",
    "\n",
    "\\begin{align}\n",
    "p_{j+1}^{n} &= exp(i(k(j+1)dx-\\omega n dt))\\nonumber\\\\\n",
    "&= exp(ik dx)\\; exp(i(kjdx-\\omega n dt))\\nonumber\\\\\n",
    "&= exp(ik dx)\\; p_j^n,\\nonumber\\\\\n",
    "\\end{align}\n",
    "\n",
    "and a grid point j-1:\n",
    "\n",
    "\\begin{align}\n",
    "p_{j-1}^{n} &= exp(i(k(j-1)dx-\\omega n dt))\\nonumber\\\\\n",
    "&= exp(-ik dx)\\; exp(i(kjdx-\\omega n dt))\\nonumber\\\\\n",
    "&= exp(-ik dx)\\; p_j^n.\\nonumber\\\\\n",
    "\\end{align}\n",
    "\n",
    "Inserting the discrete pressure wavefield solutions $p_j^{n+1}$, $p_j^{n-1}$, $p_{j+1}^{n}$ and $p_{j-1}^{n}$ in eq. (1), we get after some minor rearrangement:\n",
    "\n",
    "\\begin{equation}\n",
    "exp(-i\\omega dt)p_j^n - 2 p_j^n + exp(i\\omega dt)p_j^n = vp_j^2 \\frac{dt^2}{dx^2}\\biggl(exp(-ik dx)p_j^n - 2 p_j^n + exp(ik dx)p_j^n\\biggr).\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Assuming that $p_j^n \\ne 0$, we can divide the RHS and LHS by $p_j^n$\n",
    "\n",
    "\\begin{equation}\n",
    "exp(-i\\omega dt) - 2 + exp(i\\omega dt) = vp_j^2 \\frac{dt^2}{dx^2}\\biggl(exp(-ik dx) - 2 + exp(ik dx)\\biggr).\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "By further dividing RHS and LHS by 2, we get:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{exp(i\\omega dt) + exp(-i\\omega dt)}{2} - 1 = vp_j^2 \\frac{dt^2}{dx^2}\\biggl(\\frac{exp(ik dx) + exp(-ik dx)}{2} - 1\\biggr).\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Using the definition \n",
    "\n",
    "\\begin{equation}\n",
    "\\cos(x) = \\frac{exp(ix) + exp(-ix)}{2},\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "we can simplify this expression to:\n",
    "\n",
    "\\begin{equation}\n",
    "cos(\\omega dt) - 1 = vp_j^2 \\frac{dt^2}{dx^2}\\biggl(cos(k dx) - 1\\biggr).\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "After some further rearrangements and division of both sides by 2, leads to:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{1 - cos(\\omega dt)}{2} = vp_j^2 \\frac{dt^2}{dx^2}\\biggl(\\frac{1 - cos(k dx)}{2}\\biggr).\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "With the relation \n",
    "\n",
    "\\begin{equation}\n",
    "sin^2\\biggl(\\frac{x}{2}\\biggr) = \\frac{1-cos(x)}{2}, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "we get \n",
    "\n",
    "\\begin{equation}\n",
    "sin^2\\biggl(\\frac{\\omega dt}{2}\\biggr) = vp_j^2 \\frac{dt^2}{dx^2}\\biggl(sin^2\\biggl(\\frac{k dx}{2}\\biggr)\\biggr).\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Taking the square root of both sides finally leads to \n",
    "\n",
    "\\begin{equation}\n",
    "sin\\frac{\\omega dt}{2} = vp_j \\frac{dt}{dx}\\biggl(sin\\frac{k dx}{2}\\biggr).\n",
    "\\end{equation}\n",
    "\n",
    "This result is quite interesting. Notice, that the amplitude of the sine functions $sin(x)$ on the LHS and RHS vary between -1 and 1. However, if the factor on the RHS\n",
    "\n",
    "\\begin{equation}\n",
    "\\epsilon = vp_j \\frac{dt}{dx} \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "is larger 1 ($\\epsilon>1$), you get only imaginary solutions, while the real part is zero. Consequently, the numerical scheme becomes unstable. Therefore, the criterion\n",
    "\n",
    "\\begin{equation}\n",
    "\\epsilon = vp_j \\frac{dt}{dx} \\le 1 \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "has to be satisfied. This very important stability criterion was first described by the german-american mathematicians Richard Courant, Kurt Friedrichs and Hans Lewy in [this paper](https://gdz.sub.uni-goettingen.de/id/PPN235181684_0100?tify={%22pages%22:[36],%22view%22:%22thumbnails%22}) from 1928. The **Courant-Friedrichs-Lewy criterion** or in short **CFL-criterion**, can also be rearranged to the time step dt, assuming that we have defined a spatial grid point distance dx:\n",
    "\n",
    "\\begin{equation}\n",
    "dt  \\le \\frac{dx}{vp_j}. \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "This criterion is only correct for the FD solution of the 1D acoustic wave equation using the 3-point spatial/temporal FD operators and an explicit time-stepping scheme (eq.(1)).\n",
    "\n",
    "More generally, we can write the Courant criterion as \n",
    "\n",
    "\\begin{equation}\n",
    "dt  \\le \\frac{dx}{\\zeta vp_j}, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "where the factor $\\zeta$ depends on the used FD operator, dimension of the problem (1D, 2D, 3D) and the overall algorithm. Even though the CFL criterion strictly depends on the P-wave velocity at a specific grid point, in most cases the maximum velocity $v_{max}$ in the medium is used to estimate a constant time step $dt$ for the whole FD modelling run:\n",
    "\n",
    "\\begin{equation}\n",
    "dt  \\le \\frac{dx}{\\zeta v_{max}}, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "where $v_{max}$ is the maximum P-wave velocity in the acoustic case or the maximum S-wave velocity for the SH-problem. While the fulfillment of the CFL criterion leads to a stable simulation, it does not guarantee accurate modelling results. \n",
    "\n",
    "The CFL criterion allows us to estimate an appropriate time step dt based on the maximum velocity in the model and the spatial grid point distance. But how do we choose the spatial gridpoint distance dx?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical grid dispersion\n",
    "\n",
    "In the modelling examples at the beginning of this Jupyter notebook, we have seen that the modelled wavefield can become subject to dispersion, when choosing a too large spatial grid point distance. The result of the von Neumann analysis can also explain this behaviour. Starting from eq. (2)\n",
    "\n",
    "\\begin{equation}\n",
    "sin\\frac{\\omega dt}{2} = \\epsilon\\; sin\\frac{k dx}{2}, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "we apply the $arcsin$ to both sides\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\omega dt}{2} = arcsin\\biggl(\\epsilon\\; sin\\frac{k dx}{2} \\biggr)\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "and multiply the result by $\\frac{2}{dt}$, we get\n",
    "\n",
    "\\begin{equation}\n",
    "\\omega = \\frac{2}{dt}arcsin\\biggl(\\epsilon\\; sin\\frac{k dx}{2}\\biggr)\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Inserting this $\\omega-k$ dependence into the definition of the phase velocity \n",
    "\n",
    "\\begin{equation}\n",
    "v_{phase} = \\frac{\\omega}{k},\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "leads to \n",
    "\n",
    "\\begin{equation}\n",
    "v_{phase} = \\frac{2}{k dt}arcsin\\biggl(\\epsilon\\; sin\\frac{k dx}{2}\\biggr).\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "As you can see, the phase velocity of the numerical FD solution is a function of the wavenumber k. Therefore, it can be subject to dispersion. To investigate this problem in more detail, we rewrite the phase velocity.\n",
    "\n",
    "With the wavenumber $k=\\frac{2 \\pi}{\\lambda}$, where $\\lambda$ denotes the wavelength, we get:\n",
    "\n",
    "\\begin{equation}\n",
    "v_{phase} = \\frac{\\lambda}{\\pi dt}arcsin\\biggl(\\epsilon\\; sin\\frac{\\pi dx}{\\lambda}\\biggr).\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "From the definition of $\\epsilon = vp_0 \\frac{dt}{dx}$, we can replace $dt$ by $dt = \\frac{\\epsilon dx}{vp_0}$ in the phase velocity:\n",
    "\n",
    "\\begin{equation}\n",
    "v_{phase} = \\frac{\\lambda vp_0}{\\pi \\epsilon dx}arcsin\\biggl(\\epsilon\\; sin\\frac{\\pi dx}{\\lambda}\\biggr).\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Introducing the **number of grid points per wavelength** $N_\\lambda = \\frac{\\lambda}{dx}$, we finally get:\n",
    "\n",
    "\\begin{equation}\n",
    "v_{phase} = \\frac{N_\\lambda vp_0}{\\pi \\epsilon}arcsin\\biggl(\\epsilon\\; sin\\frac{\\pi }{N_\\lambda}\\biggr).\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Let's plot this result for $N_\\lambda$ between 2 and 12, the homogeneous P-wave velocity $vp0\\;=\\;333\\;m/s$, and $\\epsilon$ values form 0.7 to 1.0 ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nwave = np.arange(2,12,0.25) # numbers per wavelength\n",
    "vp0 = 333.0 # P-wave velocity (m/s)\n",
    "\n",
    "def dispersion_1D(eps):\n",
    "    \n",
    "    vp_phase = (vp0*Nwave/(np.pi*eps)) * np.arcsin(eps*np.sin(np.pi/Nwave)) \n",
    "    \n",
    "    return vp_phase"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vp_eps_1 = dispersion_1D(1.0)\n",
    "vp_eps_2 = dispersion_1D(0.9)\n",
    "vp_eps_3 = dispersion_1D(0.8)\n",
    "vp_eps_4 = dispersion_1D(0.7)\n",
    "\n",
    "plt.plot(Nwave, vp_eps_1, 'b-',lw=3,label=r\"$\\epsilon=1$\")\n",
    "plt.plot(Nwave, vp_eps_2, 'r-',lw=3,label=r\"$\\epsilon=0.9$\") \n",
    "plt.plot(Nwave, vp_eps_3, 'g-',lw=3,label=r\"$\\epsilon=0.8$\") \n",
    "plt.plot(Nwave, vp_eps_4, 'k-',lw=3,label=r\"$\\epsilon=0.7$\") \n",
    "plt.title('Grid dispersion')\n",
    "plt.xlabel('Number of grid points per wavelength')\n",
    "plt.ylabel('Phase velocity $v_{phase}$, m/s')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice, that no grid dispersion occurs in the case of $\\epsilon=1$. Keep in mind though that this only true for the homogeneous medium. Realistic modelling problems have a variable P-wave velocity, so we have not a constant $\\epsilon$ within the model.\n",
    "\n",
    "For all values $\\epsilon<1$, numerical dispersion can occur, if the sampling of the spatial model is too small, especially when using only the Nyquist criterion $N_\\lambda = 2$. For the 1D acoustic wave equation, the dispersion is minimized for $N_\\lambda$ values between 8-12.\n",
    "\n",
    "More generally, we can define the **grid dispersion criterion** for the spatial gridpoint distance\n",
    "\n",
    "\\begin{equation}\n",
    "dx \\le \\frac{\\lambda_{min}}{N_\\lambda} = \\frac{v_{min}}{N_\\lambda f_{max}},\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "where $N_\\lambda$ depends on the used FD operator, numerical scheme and also wave type, $v_{min}$ is the minimum P- or S-wave velocity in the model and $f_{max}$ the maximum frequency of the source wavelet. \n",
    "\n",
    "Finally, let's apply the dispersion and stability criteria to our test problem in order to find optimum dt and dx values ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# calculate dx according to the grid dispersion criterion\n",
    "Nlam = 12              # number of grid points per wavelength\n",
    "fmax = 50.0            # fmax = 2 * f0 (Hz)\n",
    "dx = vp0 / (Nlam*fmax) # spatial gridpoint distance (m)\n",
    "print('dx = ', dx)\n",
    "\n",
    "# calculate dt according to the CFL criterion\n",
    "dt = dx / vp0 # time step (s)\n",
    "\n",
    "# check CFL criterion\n",
    "epsilon = vp0 * dt / dx\n",
    "print('epsilon = ', epsilon)\n",
    "if(epsilon>1.0):\n",
    "    print('Warning: CFL condition is violated!')\n",
    "print('dt = ', dt)\n",
    "\n",
    "FD_1D_acoustic(dt,dx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What we learned:\n",
    "\n",
    "* Estimation of the Courant-Friedrichs-Lewy (CFL) stability criterion $dt  \\le \\frac{dx}{v_{max}}$ for the 1D acoustic wave equation using the von Neumann analysis\n",
    "* Dispersion analysis of the 1D acoustic wave equation\n",
    "* Grid dispersion criterion: $dx \\le \\frac{v_{min}}{N_\\lambda f_{max}}$\n",
    "* Optimize FD modelling parameters by using the grid dispersion and CFL conditions"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
